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On the Onset of Stochasticity in ACDM Cosmological Simulations 
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ABSTRACT 

The onset of stochasticity is measured in ACDM cosmological simulations using a set of clas- 
sical observables. It is quantified as the local derivative of the logarithm of the dispersion of 
a given observable (within a set of different simulations differing weakly through their initial 
realization), with respect to the cosmic growth factor. In an Eulerian framework, it is shown 
here that chaos appears at small scales, where dynamic is non-linear, while it vanishes at larger 
scales, allowing the computation of a critical transition scale corresponding to ~ 3.5Mpc//i. 
This picture is confirmed by Lagrangian measurements which show that the distribution of 
substructures within clusters is partially sensitive to initial conditions, with a critical mass 
upper bound scaling roughly like the perturbation's amplitude to the power 0.15. The cor- 
responding characteristic mass, Af cr j t = 2 10 13 M Q , is roughly of the order of the critical 
mass of non linearities at z = 1 and accounts for the decoupling induced by the dark energy 
triggered acceleration. 

The sensitivity to detailed initial conditions spills to some of the overall physical prop- 
erties of the host halo (spin and velocity dispersion tensor orientation) while other "global" 
properties are quite robust and show no chaos (mass, spin parameter, connexity and center 
of mass position). This apparent discrepancy may reflect the fact that quantities which are 
integrals over particles rapidly average out details of difference in orbits, while the other 
observables are more sensitive to the detailed environment of forming halos and reflect the 
non-linear scale coupling characterizing the environments of halos. 
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1 INTRODUCTION 

Concerns regarding the predictability of cosmological measure- 
ments in simulations have been with us for some time. In the neigh- 
bouring field of secular galactic evolutio n, it has been known for 
a while dSellwood and Wilkinson! d 19931) ) that the significant un- 
dersampling of resonances could mislead the dynamical evolution 
of N-body systems when the evolution time becomes large com- 
pared to the local dynamical time. Over the course of t he las t 
decade, various "univ e rsal" relationships dNavarro et al. I d 19971) . 
IZhang and FalU d 19991) . iRicher et al] dl99ll) ) have been extracted 
numerically from cosm ological N-body sim ulations. In this con- 
text, significant efforts dPower et alj d2003l) ) have been invested 
in comparing different numerical schemes and codes, but, with 
th e development of very high resolu t ion "z o om" re-simulations 
Weil et ail dl998l). iDiemand et alj d2004l). lHansen and Moorel 



d2006d)ISales et alj d2007h . lStrigari et alj d2007al)~ one question re- 
mains: how sensitive is a given run with respect to its initial con- 
ditions? In particular, what set of observables is likely to be robust 
with respect to a specific choice in the "phases" of the draw (the 
whitened initial realization)? In the context of cosmology, the gen- 
eral assumption has been that, even though the detailed orbits of 



dark matter particles are likely to be poorly resolved by the numer- 
ical schemes implemented, the properties of structures would nev- 
ertheless be well represented statistically. In other words, so long 
as the simulated region was large enough to represent a fair sam- 
ple of dynamically independent regions, the stochastic exponential 
departure from the unperturbed trajectories was expected to aver- 
age out when considering such a statistical sample. The question 
remains for features specific to a given realization, such as the rel- 
ative position of objects. 

The sensitivity of the gravitational N-body problem to 
small changes in initial conditions has been investigated in 
details by Kandrup and collaborators in a series of pap ers 
dKandrup and Smith! d 1 99 ll 1 1 992h ; iKandrup et ail d 1 9921 . 1 1 994h ) in 
the context of a Newtonian (time-independent) Hamiltonian. They 
have shown in particular that the growth of small perturbations in 
initial conditions is exponential, with a mean e-folding time that is 
asymptotically independent of the number of particles at large N, 
and a distribution of e-folding times that is reproducible from sim- 
ulation to simulation for sufficiently large N. In the cosmological 
context, the N-body description is an approximation of the colli- 
sionless Boltzmann equation for the evolution of the dark matter, 
so that another related question is in which sense a limit to the 
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continuum can be established as the number of particles increases . 
Indeed, it has been argued in the literature dKandrud dl99Ch : 
lGoodmanet~ai1d 19931) ) that the discretization of smooth, and pos- 
sibly integrable potentials invariably leads to strongly chaotic or- 
bits in the N-body framework, independently of the number of par- 
ticles; t his has been confirmed num erically dKandrup and Siderisl 
d200ll) : ISideris and Kandrupl d2002l) ). both for integrable and non- 
integrable underlying distributions by evolving orbits in "frozen-N 
body" samplings of the smooth mass distributions. However, Kan- 
drup and Sideris also showed that when one follows the deviation 
of orbits evolved in the frozen-N and smooth potentials with identi- 
cal initial conditions, and when the deviation amplitude is allowed 
to reach large fractions of the system size (macroscopic view), a 
continuum limit can be well defined in the sense that these macro- 
scopic departures from the orbits of the smooth potentials (which 
can be themselves either regular or chaotic) grow as a power law 
with time, and that the characteristic divergence time is growing 
with the number of particles. T hese results, together with the claim 
of universality (halo profile s l lNavarro et alj dl997|)). their shape 
Hansen and Moor el < |2006bl) ). he mass functions dZhang and Falll 



1999h , lRicher et alj ]l99l|ffl that is widely used in the cosmology 



community, lead us to revisit the problem of the sensitivity of N- 
body simulations to slight changes in the initial conditions at fixed 
power spectrum in the cosmological context, and to numerically 
investigate the presence (or absence) of "chaotic" behaviour of dif- 
ferent statistical quantities derived from N body simulations. Our 
focus will be on the transition between large scale linear dyn am- 
ics and small scale stochastic properties dStrigari et alj d2007bl) ). A 
possible concern in this context is the development of stochastic- 
ity induced by the ill-conditioning/non-linearity of the estimator of 
the chosen set of observables. Another concern lies in the specifici- 
ties of the numerical code used. Finally, numerical noise induced 
by round-off errors should also be kept at bay, since they by them- 
selves will lead to some level of stochasticity. Since the topic of 
this paper is not optimal estimation, no attempt will be made to 
argue that the set of estimators used here is superior or offers a bet- 
ter trade-off in bias ve rsus variance. Similarly, a standard integrator 
dSpringel etal.ld200ll) ) is used to carry the simulations with a set of 
conservative parameters. Round-off errors are assumed to be irrele- 
vant. Specifically, this paper will investigate what scale and mass is 
expected to play a role and will identify which quantities are found 
to be robust with respect to such exponential divergence; it will also 
find out if stochasticity breaks in as soon as non linearities occur or 
if it is possible to identify two distinct time scales in the dynamics 
of large scale structures. 

This paper is organized as follows: in Section [2] the method 
to characterize the statistical onset of stochasticity in numerical N- 
body simulations is presented. In Section[3]fhe corresponding Lya- 
punov exponents are computed for Eulerian quantities and (in Sec- 
tion [4} for Lagrangian ones. Section [5] discusses other issues and 
wraps up. 



2 METHOD & SETTINGS 

In this paper, we address the problem of the sensitivity of N-body 
simulations to initial conditions. To do so, we choose to study how 
slight changes in these initial conditions affect the evolution of the 
dispersion of a number of statistical quantities with time. This is 
achieved by generating several realizations of identical simulations 
where a small amount of random noise has been added to the initial 
realisation. The generic procedure goes as follows: 



(i) Generat e a cosmologica l simulation usin g grafic 

dBertschingen d 19951) ) and gadget dSpringel et ail d200ll) ). 

(ii) Start with the same noise file (i.e. the "phases"), but add a 
Gaussian white noise with RMS 1/30" 1 of the previous white noise 
(except in Sections |3.2| and \4.4\ where this amplitude is varied). 
This only affects the relative positions of clumps, not their spec- 
tral distribution (the expectation of the power spectrum remains 
unchanged). 

(iii) Rerun the simulation with the new white noise; 

(iv) Iterate the above two steps -50 times; 

(v) Compute a set of observables in each simulation; 

(vi) Compute the RMS (or the relative RMS) of the distribution 
of observables for various expansion factors. 

(vii) Fit the corresponding evolution of the log RMS vs the ex- 
pansion factor. 

(viii) Possibly find the scaling of its corresponding Lyapunov 
exponent (see below), with the smoothing scale associated with the 
observable (see Sec l3.2t . or the corresponding mass (see Sec l4.4t . 



Let us define the "Lyapunov exponent", Xx, as the rate of change 
of the logarithm of the fluctuation of the relevant quantity, X, as a 
function of the scale factor, a: 



A, 



din ax 
da 



(1) 



This stochasticity parameter is not strictly speaking a Lyapunov 
exponent since it corresponds neither to an asymptotic limit at 
large time, nor to an asymptotic limit at small fluctuation. It is 
closer in spiri t to th e short time Lyapunov exponent defined by 
lKandrupetai]dl997l) . 

In practice two distinct sets of simulations are considered in 
this paper, one composed of 65 realisations of 128 3 particles each 
(Si hereafter) and the other of 27 realisations with 256 3 particles 
each (1S2 hereafter). The box size is 100/i -1 Mpc, the cosmology 
a standard ACDM model (tt m = 0.3, fi A = 0.7, H = 70), 
the softening parameter is 39.5/i _1 kpc and the expansion factor 
ranges from 0.05 up to 1 for Si and from 0.05 up to 0.4 for S2 ■ 
These two sets allow us to check the robustness of our finding with 
respect to resolution. Lyapunov exponents will also be expressed 
as characteristic timescales, r, using the relationship between time 
and expansion factor in a CDM model (or equivalently in a ACDM 
model below a ^ 0.5), a oc r 2//3 . Note that the resolution in mass 
of FOF halos containing more than 100 particles corresponds here 
to 4 10 12 M G for the set Si and 5 10 11 M for the set S 2 . 



3 EULERIAN EXPONENTS 

In this Section, we investigate the "global" chaos in the evolution 
of the Eulerian properties of the density field with respect to the 
expansion factor a, as opposed to chaos in the Lagrangian proper- 
ties of objects which are specific to the matter distribution in the 
universe (such as halos and filaments). This will be addressed in 
Section|4] 



3.1 Chaos in density fluctuations 

In order to study the density fluctuations, the density fields of Si 
and S2 are sampled on a 64 3 grid using a simple NGP (nearest grid 
point) method allowing the computation of statistical quantities on 
the resulting grid, such as the average density or the density flue- 
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tuations0 Within each set, for every pixel we compute the PDF of 
the realisations of the density values at that pixel and then average 
the individual pixel PDFs over all the pixels that have mean density 
value across realisations above the given threshold. The evolution 
of the width of this pixel-averaged PDF is then computed as a mea- 
sure of the chaotic divergence amongst realisations with slightly 
different initial conditions. Specifically, Figure [T]presents the evo- 
lution of the mean relative dispersion of the density, Sp/p (where 
p is the mean density of the pixel over all the realisations and not 
the average density of the simulation), in identical pixels of the 
different realisations of Si (middle panel) and 5*2 (bottom panel), 
considering regions where density is greater than given thresholds^ 

As expected, these measurements show that this dispersion in- 
creases with time, as can easily be seen on the top panel of figure[T] 
where the PDF of Sp/p is plotted for different values of a. The fact 
that the growth rate of the dispersion increases when considering 
regions of higher densities may be explained by the higher level of 
nonlinearity of the evolution of matter distribution in these regions. 
In fact, in denser regions, the evolution becomes non-linear ear- 
lier, which favors the development of chaos. But at later times, non 
linearities have had time to develop at all considered overdensity 
levels, which explains the asymptotic merging of the curves. The 
exponential growth of the dispersions demonstrates that the evo- 
lution is chaotic as defined in the introduction and allows for the 
computation of Lyapunov exponents, Ap, as the rate of change of 
the logarithm of the average relative density fluctuation as a func- 
tion of the scale factor. 

The fact that the non-linearity in the evolution increases chaos 
is illustrated by Figure [2] where maps of the average density 
(top panel) and the corresponding Lyapunov exponent Ap (bot- 
tom panel) are plotted. Each map represents the projection of a 
10h~ Mpc slice from a sample of S2 at a = 0.35. The correlation 
between the two maps confirms the dependence of chaos on over 
density (see also the projections of different realisations of the same 
halos on Figure [6] where substructures are clearly different even 
though the shape of main halos remains mostly the same). These 
results must nonetheless be interpreted with care as the use of a 
finite sampling grid may bias the measurements. Indeed, consider- 
ing higher density regions amounts to considering smaller scale re- 
gions, of order the size of the grid pixels (« 1.5/i -1 Mpc 3 ), which 
may affect the measured value of Ap . 



3.2 Chaos transition scale 

Transition to chaotic behaviour of the density field that started with 
linear evolution is fundamentally linked to the development of the 
nonlinearity. Since different scales enter nonlinear regime at dif- 
ferent epochs, one expects that at a given time there exist a tran- 
sition scale, L c , below which variation of the density in pixels of 
the sampled field is clearly chaotic. Figure[3]presents the behaviour 
of the average value of Ap for different perturbation amplitudes A, 
as a function of the scale L. These measurements are derived by 
computing the average Lyapunov exponents in pixels on the sam- 
pled maps shown in figure [2] smoothed using a Gaussian kernel of 



1 we also considered a 128 3 grid and found no difference in the measured 
exponents. 

2 note that the number of pixels above a given threshold is going to depend 
on redshift, but for the contrasts considered here, the error on the dispersion 
due to shot noise is always negligible, as we have at least 8000 particles 
above the highest threshold, at the highest redshift. 
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Figure 1. top: Evolution of the pixel-averaged PDF of the density val- 
ues while sampling Si on a 64 3 grid for different values of a G 
[Q.l(light), ■ ■ ■ 0A(dark)] (only the pixels where p/p > 2 were con- 
sidered, where p is the cosmic mean density). As expected the full width 
half max (FWHM) of the distribution increases exponentially with the ex- 
pansion factor reflecting the chaotic behaviour of the PDF of the density 
field; middle: temporal evolution of the dispersion in the sampled density 
field per unit of the mean density, for sub regions of Si corresponding to 
thresholds in overdensity p/ p of 0.5, 1, 1.5 and 2 respectively as labelled. 
The asymptotic merging for different thresholds reflects the fact that at later 
times, regions of different overdensity levels are all in the nonlinear regime; 
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Figure 2. Logarithm of the projected density of the pixels (top frame) 
and their associated Lyapunov exponent (bottom frame), for a projected 
10h~ 1 Mpc slice S2 at a = 0.35. The comparison of the two maps empha- 
sizes the correlation of the two fields: denser regions have larger Lyapunov 
exponents. On closer inspection, one may argue that larger Lyapunov expo- 
nents lie in the outskirts of the denser regions. 

FWHM L, and considering only the overdense regions (p/p > 1). 
Density is computed by making a histogram of particles in the grid 
using the NGP method, and by smoothing it with a Gaussian kernel 
afterwards. The measurements are performed at the present time, 
a = 1 in Si simulation and at a = 0.4 for 5*2 set. 

The plot demonstrates a rather sharp transition to chaotic be- 
haviour at scales below the critical smoothing length L c ~ 3.5/i~ 
Mpc with Lyapunov exponent increasing for ever smaller scales, 
whereas on larger scales the Lyapunov exponent is small and con- 
stant. This behaviour is indicative of the ACDM background cos- 
mology of the standard model. Indeed, in the pure CDM cosmology 
with the critical density of the matter, the gravitational clustering 
would have continued to escalate to present time and one expect 
to see Lyapunov exponent falling smoothly to L ~ 8/i _1 Mpc, 
the present-day nonlinear scale Q In contrast, in ACDM cosmol- 

3 The nonlinear scale is usually defined with top-hat smoothing as 
<t 2 (-Rth) = 1. The FWHM of the Gaussian smoothing filter L that we 
use gives similar variance to the top-hat filter at i?TH ~ 0.9L. Our simula- 
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Figure 3. Evolution of the average Lyapunov exponent of the pixel den- 
sity fluctuations, Ap, as a function of the smoothing length L for regions 
where p/p > 1 and for different amplitudes of the initial perturbations, A 
(expressed as a fraction of the initial dispersion amplitude) measured in the 
set Si, The top dark dashed line corresponds to the set S2 for A = 1/30. 
The sharp transition near L S mooth ~ 3.5Mpc//i is exhibited in both res- 
olutions. The perturbation amplitude does not affect the result significantly. 
The difference in time of measurement for the S2 curve (slightly earlier 
than the freeze-out time around z ~ 1) may explain small the difference 
in the corresponding non linear scale. The bottom light dashed line corre- 
sponds also to A = 1/30 in 5i but measured on a 128 3 grid; it shows that 
the exponent is not sensitive to the sampling resolution. 

ogy the hierarchical clustering saturates when the dark energy be- 
gins to accelerate the expansion of the Universe. Numerical simula- 
tions show that in the standard ACDM model the clustering largely 
ceases by z ~ 1 ( dHatton et alj2003h ). The non-linear scale at this 
redshift is L = 3.7/i _1 Mpc, which corresponds to the mass scale 
M « 2 X 10 13 M Q . The halos of smaller mass collapse en masse 
at earlier times passing by z — 1 through a period of hierarchi- 
cal mergers with similar-mass halos as well as accretion that con- 
tributes to the formation of the chaotic features. Whereas the larger 
overdense patches, even the rare ones that turned around by z ~ 1 
and will collapse by the p resent time, evolve in a quiescent environ- 
ment of frozen hierarchy dvan den Bos ch 2002; Aubert and Pichor] 
120071) . This argues for L w 3.7/t x Mpc providing the fixed critical 
length between chaotic and regular regimes for all 2 < 1, which is 
in general agreement with our measurements. 



4 LAGRANGIAN EXPONENTS 

In the previous section, we studied the development of chaos in the 
density field of cosmological simulations. We measured the evolu- 
tion of the variance of this density field on a grid (i.e. at peculiar 
Eulerian locations) and showed that chaos tends to be more pro- 
nounced in higher density regions as well as on smaller scales. Let 
us now focus instead on Lagrangian properties of peculiar objects 
with a physical significance such as dark matter halos or filaments. 

tions are normalized to ff(8fc _1 Mpc) = 0.92 which at a = 1 corresponds 
to nonlinear scale Rth ~ 7.2h~ 1 Mpc, i.e L R3 8/i _1 Mpc. 
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4.1 Inter skeleton distance 

Filaments correspond to a central feature of the large scale distri- 
bution of matter: large void regions are surrounded by a filamen- 
tary web linking haloes together. Studying the properties of the 
filaments isn't an easy thing and one first needs to find a way of 
extracting their location from a simulation. The skeleton gives a 
mathematical definition of the filaments as the locus where, start- 
ing from the filament type saddle points (i.e. those where only one 
eigenvalue of the Hessian is positive), one reaches a local maxi- 
mum of the field by following the gradient. This is equivalent to 
solving the equation: 



for x, where p(x) is the density field, V/9 its gradient, and x the 
position. Although apparently simple, solving this equation is quite 
difficult which is why a lo cal approximation was introduced in 
dSousbie et all < |2006|. 120071) ): the local skeleton. One can show 
that, up to a second order approximation, solving Equation {2]l is 
equivalent to finding the points in the field where the gradient is 
an eigenvector of the Hessian matrix together with a constraint on 
the sign of its eigenvalues. This approach leads to a system of two 
differential equations, solved by finding the intersection of two iso- 
density surfaces of some function of the density field and its first 
and second derivatives. This procedure is very robust and allows 
for a fair detection of the dark matter filaments. Figure [4] displays 
the skeleton of different realisations of 5*2 at a(t) = 0.1 (top) and 
a(t) = 0.4 (bottom). Note that the dispersion of the skeleton lo- 
cation has increased with the scale f actor. Using the method de- 
scribed in JSousbie et aiH200dl2007l) ). the local skeleton provides 
a list of small segments. In order to measure the distance between 
two skeletons, for each segment, the distance to the closest seg- 
ment in the other skeleton is computed leading to the PDF of this 
distribution. The mean distance between the two skeletons is set 
to t he position of the fir st mode of their inter-distance PDF (See 
also ICaucci et all d2008h ). We then define a mean inter- skeleton 
distance among all the realisations within a set as the arithmetic 
average of their pairwise distances. This means that the normal- 
ized inter-skeleton distance, (D) /Lo, is our measure of the dis- 
persion in the skeleton location. It is a Lagrangian property since 
it follows the flow. Its evolution as a function of the scale factor 
is plotted on Figure [5] for different smoothing lengths Lo, for £2 
(top) and Si (bottom). The smoothing operation is achieved, as 
previously, by convolving the density field with a Gaussian ker- 
nel of FWHM Lo, ranging from Lo = 1.2 /i _1 Mpc (3 pixels) up to 
Lo = 3.5 h~ 1 Mpc (9 pixels). It is clear that whatever the smooth- 
ing scale or the resolution used, the evolution of the dispersion is 
linear with the scale factor. A shift in the skeleton of the initial 
conditions will evolve linearly with time and not exponentially: the 
skeleton at present time won't be affected very much. There is no 
chaotic drift of the position of the skeleton and thus no chaos in 
the evolution of the cosmic web. Note nonetheless that the smaller 
the smoothing length, the stronger the increase of (D) / Lo. This 
implies that smaller scales are more sensitive to initials conditions, 
which is confirmed by the fact that (D) /Lo is larger for lower val- 
ues of Lo, whatever the value of a. 



4.2 Positions of halos 




X (h- 1 Mpc) 




Figure 4. The local skeletons of the different realisations of ^2 , computed 
at a(i) = 0.1 (top) and a(t) = 0.4 (bottom) and for a smoothing length 
Lo = 1.2/i —1 Mpc. Each figure corresponds to the projection of a 10h~ 1 
Mpc thick slice. Each color represents a different realisation of the simula- 
tion, the color coding is not consistant between the top and the bottom pan- 
els. The dispersion in the position of the skeletons appears to have grown 
from a(t) = 0.1 to a(t) = 0.4. 



Turning to stochasticity on smaller scales in a Lagrangian frame- 
work (i.e. ignoring the absolute shift in position relative to a fixed 
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Figure 5. The normalized mean distance, (D) /Lq, between the skeletons 
of S2 (top) and Si (bottom) as a function of the scale factor and for different 
values of the smoothing length, Lq in h~ 1 Mpc. Because of the lack of 
accuracy at smaller scales, only the larger smoothing lengths are represented 
for Si. At these resolutions, the two sets agree. The cosmic web dispersion 
clearly evolves linearly with time, confirming that chaos is linked to non- 
linearities. 
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Figure 6. Heaviest cluster (in the X-Y plane in Mpc/h) of 9 realisations 
of S2 at z = 1.5. The position and the global shape of the halo does 
not change from one simulation to another, but the substructures are quite 
different; this is confirmed via automated substructure identification using 
ADAPTAHOP. 



frame), let us define a matching procedure to identify structures in 
different runs, halo e s are first identified using th e FOF algorithm 
dDavis et al. I dl985h ; ISugInohara and Sutol Jl992h ) with a percola- 
tion length of 0.25/i _1 Mpc for Si and 0.5/i _1 Mpc for 5*2 corre- 
sponding to 0.2 x mean interparticular distance. In order to tag dif- 
ferent FOF haloes in different realisations as counterparts, all par- 
ticles of a given halo are matched in another realisation using their 
initial index (Figure[6]l. The halo of the other simulation containing 
most of these particles is tagged as its counterpart. The procedure is 
carried over all pairs of simulations, allowing the measurement of 
the variation in the halo properties like their spins, their positions, 
their velocity dispersion tensors or their masses. 

As shown on Figure [6] the haloes locations seem relatively 
insensitive to small changes in the initial conditions. The evolu- 
tion of the mean distance between a halo in a given simulation and 
the same halo in another realisation is linear, as for the skeleton, 
which confirms the first impressions: no chaos is observed at linear 
scales and so Xq, the Lyapunov exponent of the inter halo distance, 
is null. But the most interesting results involve the substructures. 
The halo pictured in Figure[6]is a good example of the generic be- 
haviour. The number of substructures changes from one realisation 
to another (here, 1 or 2 substructure(s)) and their positions also dif- 
fer. These results are confirmed by an automated detec tion of the 
substructures using ADAPTAHOP l lAubert et all 1 120041) ). Both the 
locations and the number of substructures are possibly subject to 
chaos, but the lack of a cross identification procedure makes it dif- 
ficult to quantify it and is somewhat beyond the resolution of these 
sets of simulations (Section[43]addresses this problem for the FOF 
halos). This trend confirms quantitatively the findings of section l3~Tl 
from the point of view of Eulerian estimators which are sensitive 
to the detailed extension of the distribution of matter within halos: 
denser regions were found to be chaotic, and will be addressed in 
more details in Section 14.41 in terms of halo density and velocity 
moments. 



4.3 Connexity and mass of clusters 

The connexity of haloes can be defined as follows: considering 
the p*' 1 halo, Hp, in the r th realisation, its particles are spanned 
amongst n haloes in the r' th realisation, and a fraction (k G 
1, .., n) of them belong to a halo k among n in the r' th realisation. 
Hence, its relative connexity Cp T can be defined as: 



(3) 



where by construction C£ r is equal to one if both haloes are iden- 
tical in realisations r and r'\ Cp T is equal to n if the halo p splits 
into n haloes with equal fractions fH. — 1/n in realisation r', 
while preserving continuity when the values of /T£ differ, see Ap- 
pendix A. The mean connexity, C, is obtained by averaging over all 
haloes containing more than 100 particles in every possible com- 
binations of realisations and is a measure of the dispersion of the 
particles. 

As shown on Figure|7] C increases with the scale factor, rang- 
ing from 1.17 (i.e., statistically, 90% of the particles belong to a 
unique halo in other realisations) to 1.37 ( 85%) from a — 0.2 up 
to a = 1.0. The connexity clearly does not vary exponentially with 
the scale factor : there are statistically no halo fission during evo- 
lution (thanks to the efficiency of dynamical friction). Moreover, 
two haloes marginally linked by FOF would almost always end up 
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Figure 7. The haloes average connexity computed over all the realisations 
of Si as a function of the scale factor. While the connexity is not subject to 
chaos, its value increases with time. This result can be understood through 
the difference in merging time of the haloes. 



merging sooner or later. More and more haloes merge at different 
times in different realisations which is in part due to the fact that 
some threshold is involved in the FOF algorithm: a precise linking 
length has to be chosen, inducing the possibility that small changes 
in particles position can induce significant changes in halo merg- 
ing time (according to the FOF definition of a halo). At later times 
(a > 0.7), the connexity reaches a plateau, which suggests that 
when haloes are massive enough (M ^ M c , see Section |4j4] be- 
low), they become insensitive to the merging of lighter ones, since 
equal mass merging rarely occur below z = 1. The analysis of the 
masses of the haloes shows that there is no sweeping change and so, 
no obvious evolution of the haloes mass distribution: the associated 
Lyapunov exponent, Am, is null. The number of different particles 
increases with time but the missing particles are replaced by new 
particles. Thus, the mass stays quite constant even as the connex- 
ity increases. It follows that the mass function extracted from the 
N-body simulations are found to be quite robust with respect to 
changes in the initial conditions. Although the masses of the haloes 
are similar in different realisations, some of the particles which 
compose them may be different, which may generate differences 
in the physical properties of the haloes. The substructures are dif- 
ferent (Fig.[6]l in their numbers and positions, which is responsible 
for the Eulerian chaos found in Section [3~2l Let us now re-explore 
this in a Lagrangian framework. 



4.4 Spin Orientation of clusters 

The influence of chaos on the spin of haloes is estimated by com- 
puting the cosine of the angle, 8 pq , between the spins J p and J q of 
corresponding haloes in two different realisations p and q: 



a — arccos 



1 c 



COS 9r> 



(5) 



cos (6 P 



J p ' J q 



(4) 



where the sum is over all the N c possible pairwise combinations of 
realisations. Note that only bins of masses containing more than 30 
haloes have been retained. 

Figure [8] displays the exponential growth of this dispersion 
with time, and shows that the precise value of its associated Lya- 
punov exponent A CT depends on the selected bin of mass. It also 
shows that the exponent does not seem to be sensitive to shot noise, 
as its value is left unchanged when resolution is increased between 
5i and S2. 

Also, as it is seen in Fig(6] the detailed distribution of satellites 
within a given cluster varies from one realisation to another; the an- 
gular momentum orientation (in contrast to say, its modulus or the 
halo mass) is quite sensitive to the outer region of the distribution. 
Recall that the spin p a rameter (i.e. t he A = J/(\/2MV2oo-R20o) 
jBullock et all J200lh . lAubert et al.l fe004h ) of a halo displays no 
chaotic behaviour. It stays quite constant from a simulation to an- 
other and the evolution of its dispersion is not exponential. 

The measured Lyapunov exponent ranges from to 0.3. The 
value of the mean dispersion of the orientation of the spin for heav- 
ier haloes is about 35 degrees (= exp(3.55)). Globally this sug- 
gests that the orientation of the spin varies with the tidal field, 
which in turn depends on the relative position of structures within 
the environment of the halo. 

For lighter haloes, the measured value of X a is higher than 
for heavier ones (Figure [8} which may be partly explained by the 
fact that a slight change in a few clumps within the haloes has a 
larger influence on its spin wh en they represent a signifi cant frac- 
tion of it. Faltenbacher & al. dFaltenbacher et all d2005h ) showed 
that if the lightest halo has a mass less than 10% of the mass of the 
larger halo, the orientation of the resulting post merging halo will 
remain statistically the same. In contrast, if its mass is greater than 
20% of the more massive halo, the final orientation of the merged 
halo depends on the speed vector of the two progenitors. Our re- 
sults corroborate well their finding since the lightest haloes that 
are formed by merging of two substructures of comparable masses 
have chaotic spins (substructures being chaotic, see Section [3] and 
Figure [6), while heavier ones have spin that are relatively stable 
with time (they only merge with much smaller haloes). It emerges 
from these measurements that there is a critical mass, M c , above 
which chaotic behaviour disappears. Haloes heavier than this mass 
are too heavy to feel the influence of incoming clumps and their 
spins are clearly defined. They are therefore not subject to chaos 
and their Lyapunov exponents are null at the one sigma level, in 
contrast to lighter ones whose spin are sensitive to the initial con- 
ditions and whose Lyapunov exponents are positive. 

As for the critical smoothing length (see Section [3), we can 
study the evolution of this critical mass as a function of the ampli- 
tude, A, of the perturbations. Figure|9]shows this evolution. A good 
fit of this transition mass is given by M c = 2 10 13 M G A - 15 . The 
higher the amplitude of the perturbations, the higher the required 
time for haloes to have a spin clearly defined. Consequently, the 
critical mass increases with the perturbation amplitude, and haloes 
that can be considered stable are heavier. Note that it means that 



4 the estimator of the dispersion, Eq. (3) is robust since weighting the sum 
by the spin parameter yields the same results. 
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Figure 8. Logarithm of the dispersion of the angle between the spin 
of one halo of Si as a function of the scale factor. Results are com- 
puted for different ranges of masses, 5 10 12 M o <M<6 10 12 M q (top) and 
1.6 10 13 M©<M<2 1O 13 M0 (bottom). For heavier halos the Lyapunov ex- 
ponent vanishes. The triangles correspond to the first class of lighter mass, 
but measured in 52 ; the exponent remains unchanged which suggests that 
particle shot noise is not an issue. 

the spin is constant with time but not very reliable since its final 
orientation depends, in part, on the initial conditions. 

4.5 Orientation of the velocity dispersion tensor 

The orientation of the velocity dispersion tensor is also a quantity 
of interest from the point of view of stochasticity since it is related 
to the shape of the halo via the Virial theorem. The correspond- 
ing estimator involves computing the orientation of the eigenvector 
V associated to the largest eigenvalue of the velocity dispersion 
tensor. As for the orientation of the spin, (Sec. I4,4t the angle 8 vq 
between the eigenvector V v of the halo in simulation p and its cor- 
responding eigenvector V q in a simulation q is computed as: 



cos(# p «j) = 



VyV q 

\v P \\\\v q \ 



(6) 



For every bin of mass, a measure of the dispersion, a, is also given 
by Equation l[5}- Once again, only bins of masses containing more 
than 30 haloes were considered. As shown in Figure \10\ this esti- 
mate is consistent with the exponents of the orientation of the spin: 
only the lightest masses are sensitive to the initial conditions, while 
the dispersion of the orientation for the heavier masses is constant 
(about 40 degrees). The measured Lyapunov exponent, At, ranges 
from up to 0.65 . These results corroborate well those for the spin 
axis, given that its orientation follows the third eigenvector of the 
disp ersion matrix (i.e. the axis along which dispersion is the small- 
est). iFaltenbacher et al] d2005T) showed that the orientation of the 
principal axis of the halo is correlated with the vector linking the 
two mergers (i.e. their relative positions) particularly in the case 
where one of the mergers has a mass smaller than 10% of the sec- 
ond one. The chaos found at small scales (substructures scale) is 
once again responsible for the chaos in these geometrical properties 
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Figure 9. Critical mass, M c , (in units of 10 10 Mq) for the spin orientation 
as a function of the amplitude of the perturbations, A (in fraction of the 
initial dispersion amplitude). It appears that M c = 2 10 13 Mq A 015 . The 
larger the amplitude of the perturbation, the heavier the haloes that can be 
considered stable. 



of haloes since changing the initial conditions amounts to changing 
the relative positions of the substructures (see sec l3.lt and thus to 
changing the orientation of the resulting halo's dispersion tensor. 
As for the spin, there is evidence of a critical mass above which 
chaotic evolution disappears: more massive haloes only merge with 
lighter ones that do not affect their global properties. 



5 CONCLUSION AND DISCUSSION 

Let us first emphasize here again that the term chaos is used in this 
paper in the loose sense, as the age of the universe does not allow 
for many e-foldings on larger scales. Table Q] summarizes the dif- 
ferent Lyapunov exponents computed in this paper. As shown in 
section |3T| (Figure [3l, chaos appears below a critical scale which 
corresponds roughly to cluster scales. The higher the density, the 
more chaotic is the corresponding region. We also found that both 
Lagrangian and Eulerian measurements are consistent: super clus- 
ters and filaments, whose dynamics is globally linear (large scale 
structures), are not stochastic: a shift in the initial conditions will 
increase linearly with the scale factor. By contrast, the distributions 
of substructures within clusters, whose characteristic size is smaller 
than ~ 3.5/i _1 Mpc, are governed by non-linear dynamics and may 
undergo a stochastic evolution for some observables. 

Nevertheless, this chaos at substructures scale does not occur 
for all physical characteristics of the cluster's halo. A main fraction 
of particles remains in the same halo from one realisation to an- 
other, while a difference arises (in part) from the delay in merging 
times of the substructures. These timing effects are however aver- 
aged out yielding, to first order, a constant halo mass. It follows 
that the mass function derived from a simulation is quite consistent 
from one realisation to another. Similarly, the dispersion of the am- 
plitude of the total spin of haloes does not increase exponentially 
with time. 

The mass of a given halo is an integrated quantity which does 
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Figure 10. Logarithm of the dispersion of the angle between the first eigen- 
vector of the velocity dispersion tensor of the haloes of Si , as a function 
of the scale factor. Results are computed for different ranges of masses: 
5 1O 12 M <M<6 1O 12 M (top) and 2.3 10 13 M Q <M<2.8 1O 13 M (bot- 
tom). The less massive haloes are more sensitive to the initial conditions, 
the average angle being constant for the heavier ones, about a value of 40 
degrees. The triangles correspond again to the 52 set and shows not differ- 
ence. The Lyapunov exponent, At ranges from to 0.65. 



not trace which specific particle entered the FOF halo; similarly, 
the spin parameter is also an adiabatic invariant, and the trace of the 
dispersion tensor will relax rapidly to its virial expectation (which 
is mass dependent) in a few short dynamical times; in contrast, the 
spin orientation or the orientation of the dispersion tensor will de- 
pend precisely on the orientation of velocities of the entering parti- 
cles and has no direct relation to the mass of the halo; it also reflects 
the initial environment of the p roto halo. Fo r instan ce it has been 
shown in lSousbie et alj ilOOH) , IXubert et al.l f2004h that the halos 
preferentially anti-align their spin with the axis of the filament in 
which they are embedded, while we have shown in Section PTTI that 
the filament's locus was not stochastic. 

It is possi ble to recast these interp retations in the context of 
the peak-patch teond and Mversl J 19961) ) description of haloes. In 
this framework, massive haloes correspond to large quasi spherical 
patches around density peaks, which non-linear evolution will de- 
couple from their neighbouring large patches thanks to the cosmic 
acceleration below z ~ 1. Conversely, small haloes correspond 
to small typically aspherical peak-patches, and will acquire tidal 
torques early on which depend specifically on the detailed white 
noise realisation (which fixes the shape of the peak-patches). In the 
tidal torque theory, the mass and the spin parameter are essentially 
integral functions over the volume of these patches, hence will not 
depend on the initial perturbations, whereas the spin orientation it- 
self is sensitive to these perturbations, at least at the lower end of the 
mass spectrum. This is consistent with th e low scatter relation ship 
between the spin parameter and the mass jAubert et af]j2004l) ). 

Thanks to angular momentum leverage, the orientation of 
haloes is itself affected by stochasticity mostly at small scales, a 
result which seems insensitive to shot noise as the lyuponov ex- 
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Table 1. Lyapunov exponents of the different observables studied. Inter- 
estingly, many global properties of halos do not display chaotic behaviour. 



ponents are consistent between sets Si and & ■ In fact, as long as 
the haloes merging together have similar sizes (masses), the ori- 
entation of both spin and velocity dispersion tensors is determined 
by the relative positions and velocities of the two mergers, whose 
dynamics is non-linear and whose characteristic size is below the 
critical scale (Figure[3)- These results seem robust with respect to 
resolution. 

When the halo is formed and well-isolated by cosmic accelera- 
tion, it merges only with satellites/substructures whose masses rep- 
resent a small fraction of the host's mass. Consequently, their ori- 
entations are globally preserved after merging, and thus, the chaotic 
behaviour stops and the dispersion in the orientation remains at the 
same level (i.e. the resulting average angle is unchanged). Hence, 
a critical mass can be defined as the mass above which this chaotic 
behaviour of the orientation stops. The measured value of this crit- 
ical mass, M c = 2 1O 13 M A 015 , is just below the scale of non- 
linearity at z ~ 1 and shows weak dependence on the amplitude 
of the added perturbative noise: M c = 2 lO 13 M yl 015 . Although 
some slow increase of M c with A is expected since adding power 
to inhomogeneities shifts the nonlinear scale to higher masses, the 
details of the dependence require further investigation. 

This paper has concentrated on a realistic ACDM cosmology: 
it would also be interesting to rerun this investigation on scale-free 
power spectra to confirm that the dark energy is indeed responsible 
for the saturation of M c . A natural extension of this work, clearly 
beyond its current scope, would also involve computing Lyapunov 
exponen ts for the properties of substructures within halos (see for 
instance rValluri et all d2007l) ). and parameters corresponding to the 
inner structure of halos, such as NF W concentrat i on par ameter, the 



phase space density Q = pp/g 3 1 Peirani et al.l | |2006|) ). the Gini 



index or the asymmetry dConselice et al.l d2007t) ) within the FOF. 

In closing, the answer to our riddle is that chaos and non- 
linearities are very strongly linked, and both occur at small scales 
(substructures scales) though some non linear halo parameters 
(spin, mass etc..) do not seem to be subject to chaos. While the 
large scale structures in a simulation (filaments and haloes) are 
quite robust both in their locus and properties, the distribution of 
substructures is more sensitive to initial conditions since their num- 
bers and positions vary when initial conditions vary. This in turn 
may prove to be a concern when generating zoomed resimulations. 
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APPENDIX A: CONNEXITY 

Let us consider a halo, H, split into n parts, P™, i ^ n, with a 
fraction, fi, of its particles in each of them, the indices i being 
sorted following the decreasing values of /; (fi fj if i < j). 
A measure C n of the connexity of H should indicate the number 
of clumps into which it was split, this number does not necessarily 
have to be an integer, depending on the fraction of the mass of H 
that went into each P™. For instance, if n = 2, we want to obtain 
C2 = 2 when fi — fi = 1/2 and Ci — > 1 when /i — » 1 and 
fi — 1 — /i — » 0; as the indices are sorted, < fa < 1/2. So, in 
this case, we could write the connexity of H as: 

C 2 = 1 + 2/2. (Al) 

Now, considering that // was split into 3 parts Pf, then < 
(/ 2 + fa) < 2/3 and < fa < 1/3. So = 2 + 3fa -> 2 
if / 3 -> and C3 -> 3 when fa -> 1/3. It follows that C3 = 
{h + /3)C3 2 when / 2 1/3 (which implies that / 3 -> 1/3) 
and that C'J — > when — » (which implies that /3 — > also). 
So 

C 3 = l + (/2+/3)(2 + 3/ 3 ) (A2) 

has the right properties to represent the connexity of a halo split 
into 3 parts. Hence, by generalizing recursively this formula, we 
obtain: 

C n = 1 + (/ 2 + •••/„) [2 + (/ 3 + •• • /„) [• • • [n - 2+ (A3) 
(/n-i + /«) [(n - 1) + nf n ]} ■■■]], (A4) 
which can be developed as: 

C n = 1 + 2(/ a + • • • jf n ) + 3(/ a + • ■ • /n)(/ 3 + •••/„)+ (A5) 
•• • + n(/ 3 + • •• + /n)(/ 3 +••• + /„)••■ (/») , (A6) 

ri fin \ 

= E 1 IIEa . ^ 

»=i y=i fc=j y 

which corresponds to Eq. ((3). Note that by construction, the braket 
in Eq. \A7\ is smaller than 1/i, so that C n is always smaller or 
equal to n. 



